#modified version of fixedSubs to work out transition/tranversion ratio for each gene
fixedSubs <- function(chick.transcript,grouse.transcript,dir="/Paterson/Datafiles/grouse/test_files",transition.ratio=FALSE){
	#calculates fixed differences between a pair of sequences
	#(called chicken and grouse)
	require(GeneR)
	muscle.char <- character(4)
	muscle.char[1] <- paste(">ChickenProtein")
	muscle.char[2] <- strTranslate(chick.transcript)
	muscle.char[3] <- paste(">GrouseProtein")
	muscle.char[4] <- strTranslate(grouse.transcript)
	
	muscle.file <- paste(dir,"/proteins_tmpfile.fa",sep="")
	muscle.file.out <- paste(dir,"/proteins_tmpfile.afa",sep="")
	writeLines(muscle.char,muscle.file)
	
	ta.char <- character(4)
	ta.char[1] <- paste(">ChickenTranscript")
	ta.char[2] <- chick.transcript
	ta.char[3] <- paste(">GrouseTranscript")
	ta.char[4] <- grouse.transcript
	
	tranalign.file <- paste(dir,"/transcripts_tmpfile.fa",sep="")
	tranalign.file.out <- paste(dir,"/transcripts_tmpfile.afa",sep="")
	writeLines(ta.char,tranalign.file)

	#need to set R up use same PATH as X11
	system(paste("/usr/local/bioinf/muscle/muscle -in",muscle.file,"-out",muscle.file.out))
	system(paste("/opt/local/bin/tranalign",tranalign.file,muscle.file.out,tranalign.file.out))
	require(Biostrings)
	prot.afa <- read.BStringSet(muscle.file.out)
	dna.afa <- read.BStringSet(tranalign.file.out)
	
	tst <- strsplit(as.character(prot.afa),"")
	rm.gaps <- unique(c(grep("-",tst[[1]]),grep("-",tst[[2]])))
	tst2 <- list(chick.prot = toupper(tst[[1]][-rm.gaps]), grouse.prot = toupper(tst[[2]][-rm.gaps]))
	
	#insert homology criteria here
	# .. nothing yet ..
	
	#yn00 estimates
	paml.estimates <- try(pamlCalc(dna.afa,dir))
	#note, if paml fails may still take results from existing files
	
	prot.diff <- sum(tst2$chick.prot!=tst2$grouse.prot)
	tst <- strsplit(as.character(dna.afa),"")
	rm.gaps <- unique(c(grep("-",tst[[1]]),grep("-",tst[[2]])))
	tst2 <- list(chick.dna = toupper(tst[[1]][-rm.gaps]), grouse.dna = toupper(tst[[2]][-rm.gaps]))
	if(any(!tst2$chick.dna %in% c("A","G","C","T"))) warning("unexpected base in chicken sequence\n")
	if(any(!tst2$grouse.dna %in% c("A","G","C","T"))) warning("unexpected base in grouse sequence\n")
	dna.diff <- sum(tst2$chick.dna!=tst2$grouse.dna)
	if(transition.ratio==TRUE){
		#count up number of transitions and tranversions
		tst2$chick.type <- sapply(tst2$chick.dna,function(X){
			if(X %in% c("A","G")) return("R") #purine
			if(X %in% c("C","T")) return("Y") #pyrimidine
			warning('unknown base in chicken sequence')
			return("N") #only if function gets this far
			})
		tst2$grouse.type <- sapply(tst2$grouse.dna,function(X){
			if(X %in% c("A","G")) return("R") #purine
			if(X %in% c("C","T")) return("Y") #pyrimidine
			warning('unknown base in grouse sequence')
			return("N") #only if function gets this far
			})
		no.transversions <- sum(tst2$chick.type!=tst2$grouse.type)
		no.transitions <- dna.diff - no.transversions
		tmp <- c(prot.diff,dna.diff-prot.diff,no.transitions,no.transversions)
		names(tmp) <- c("Dn","Ds","transitions","transversions")
		}else{
		tmp <- c(prot.diff,dna.diff-prot.diff)
		names(tmp) <- c("Dn","Ds")
		}
	if(class(paml.estimates)=="try-error"){
		tmp <- c(tmp,rep(-1,9))
		}else{
		tmp <- c(tmp,unlist(paml.estimates))
		}
	tmp
	}


#for(grouse.gene in c(18))){#change to nrow(grouse.chick.map)
grouse.gene <- 9
	chick.transcript <- as.character(chick.seqs[[grep(grouse.chick.map$ensembl_transcript_id[grouse.gene],names(chick.seqs))]])
	grouse.transcript <- getCDS(seq.name=grouse.chick.map$grouse.id[grouse.gene],single.cds=TRUE)
	if(!"transcript" %in% names(grouse.transcript)) next
	
	#remove *s from protein/transcript
	star.pos <- grep('\\*',strsplit(grouse.transcript$protein,'')[[1]])
	if(length(star.pos)>0){
		tst.trans <- strsplit(grouse.transcript$transcript[[1]],'')[[1]]
		cut.trans <- tst.trans[-c( (star.pos-1)*3+1, (star.pos-1)*3+2, (star.pos-1)*3+3 )]
		grouse.transcript$transcript <- paste(cut.trans,collapse="")
		}
	
	test.MK <- try(fixedSubs(chick.transcript,grouse.transcript$transcript,transition.ratio=TRUE))
	if(class(test.MK)=="try-error") next
	grouse.chick.MK[grouse.gene,c("Dn","Ds","transitions","transversions")] <- test.MK
	}

#useful code to set PATH variable
Sys.setenv(PATH=paste(Sys.getenv("PATH"),"/usr/local/bioinf/paml/bin",sep=":"))

#use PAML to get transition/transversion for each sequence plus dn/ds?
pamlCalc <- function(sequences,dir){
	require(Biostrings)
	#sequences is a DNAstringset object
	#write sequences into phylip sequential format
	paml.file <- paste(dir,"/paml_input.phy",sep="")
	if(length(sequences[[1]])!=length(sequences[[2]])) stop('sequences of unequal length\n')
	paml.char <- character(length=1+2*length(sequences))
	paml.char[1] <- paste(length(sequences),length(sequences[[1]]))
	for(i in 1:length(sequences)){
		paml.char[2*i] <- names(sequences)[i]
		paml.char[1+2*i] <- toupper(as.character(sequences[[i]]))
		}
	writeLines(paml.char,paml.file)
	#assumes that there are control files yn00.ctl and codeml.ctl in dir
	#outputs to paml_output.out
	dir.now <- getwd()
	setwd(dir)
	system("yn00") #run Yang and Neilsen 200
	paml.grep <- system("grep -A 2 \'kappa\' paml_output.out",intern=TRUE)[c(1,3)]
	#1     2       3       4        5     6       7        8     9    10    11
	#seq. seq.     S       N        t   kappa   omega     dN +- SE    dS +- SE
	paml.out <- as.numeric(strsplit(paml.grep[2]," ")[[1]][grep('^[0-9]',strsplit(paml.grep[2]," ")[[1]])])
	setwd(dir.now)
	list(dN=paml.out[8],dN.se=paml.out[9],dS=paml.out[10],dS.se=paml.out[11],
		kappa=paml.out[6],omega=paml.out[7],Ls=paml.out[3],Ln=paml.out[4],t=paml.out[5])
	
	
	}